\section{Code for problem 4}
\subsection{Code for problem 4.a}
\label{code4a}
\small
\begin{verbatim}
# This function generates a severity of an earthquake.
generateEarthquakeSeverity <- function()
{
  # This is the inverse distribution function for the jump sizes.
  inverseDitribution <- function(u)
  {
    return(2/sqrt(1-u) - 1)
  }

  # The threshold is set at 20.
  s <- 20
  S <- 0

  # We continue as long as we haven't exceeded the threshold.
  while(S < s)
  {
    # We generate another jump in the compound Poisson process.
    u <- runif(1)
    S <- S + inverseDitribution(u)
  }

  # We return the amount with which we exceeded the threshold.
  return(S-s)
}

# The size of the sample that we want to generate.
N <- 1000

# Generating the sample.
sample <- c()
for(i in 1:N)
{
  sample <- c(sample, generateEarthquakeSeverity())
}

# Drawing the histogram of the generated sample.
hist(sample, breaks = 200)
\end{verbatim}

\subsection{Code for problem 4.b}
\label{code4b}
\small
\begin{verbatim}
simulateJ2 <- function(time)
{
  rootFunction <- function(y, a)
  {
    1 - (y^(-1/4)+y^(-3/4))/2 - a
  }

  t <- 0
  j_2 <- 0
  while(t < time)
  {
    t_i <- rexp(1, 8/5)
    t <- t + t_i
    if(t < T)
    {
      u <- runif(1)
      xi_i <- uniroot(rootFunction, c(1, 1000000), a = u)$root
      j_2 <- j_2 + xi_i
    }
  }
  
  return(j_2)
}

simulateJ1 <- function(time)
{
  x <- rnorm(1, 0, sqrt(time))
  return((x + (3/10)*time)*sqrt(3/5))
}

r <- 0.04
sigma <- 0.3
S_0 <- 100
T <- 1
K <- 105

# The size of the sample that we want to generate.
N <- 10000

mu <- 0
for(i in 1:N)
{
  w <- rnorm(1, 0, sqrt(T))
  j1 <- simulateJ1(T) 
  j2 <- simulateJ2(T) 

  S_T <- S_0*exp((r-(sigma^2)/2)*T + sigma*w + j1 + j2)
  
  mu <- mu + max(S_T - K, 0)
}

mu <- mu*exp(-r*T)/N

mu
\end{verbatim}

\newpage
